#版本1####
library(monocle)
# BiocManager::install("monocle")#,force = TRUE
setwd("C:/Users/ZFB/Desktop/单细胞生信/GSE212966")
Only_T<-readRDS('./data/temp/Only_T.rds')
Only_T_matrix<- as(as.matrix(Only_T@assays$RNA@counts),'sparseMatrix')
Only_T_data<- Only_T@meta.data
Only_T_data$celltype<-Only_T@active.ident
Only_T_f_data <- data.frame(gene_short_name = row.names(Only_T),row.names = row.names(Only_T))
pd<-new('AnnotatedDataFrame',data=Only_T_data)
fd<- new('AnnotatedDataFrame',data=Only_T_f_data)
#你可以去看一下identical(rownames(fd),rownames(Only_T_matrix))返回的结果是否为TURE，若不是可以运行下面的代码
Only_T_matrix<- Only_T_matrix[rownames(fd), ]

cds<- newCellDataSet(Only_T_matrix,
                     phenoData = pd,
                     featureData = fd,
                     lowerDetectionLimit =0.5,
                     expressionFamily =negbinomial.size())
cds<-estimateSizeFactors(cds)
cds<-estimateDispersions(cds)
cds <- detectGenes(cds,min_expr = 0.1)
print(head(fData(cds)))
expressed_genes <- row.names(subset(fData(cds),num_cells_expressed >=10))  #过滤掉在小于10个细胞中表达的基因，还剩11095个基因:

express_genes<-VariableFeatures(Only_T)
cds<-setOrderingFilter(cds,express_genes)

p1 <- plot_ordering_genes(cds)
p1

diff<-differentialGeneTest(cds[expressed_genes,],fullModelFormulaStr = "~celltype",cores = 1)
deg<-subset(diff,qval<0.01)
deg<-deg[order(deg$qval,decreasing = F),]
write.table(deg, file = "./data/temp/monocle2_deg.xls",col.names = T,row.names = F,sep = "\t",quote = F)

ordergene<-rownames(deg)
cds<-setOrderingFilter(cds,ordergene)

p2 <- plot_ordering_genes(cds)

cds<-reduceDimension(cds,max_components = 2, method='DDRTree')
cds<- orderCells(cds) # In dfs(graph = graph, root = root, mode = mode, unreachable = unreachable,  :Argument `neimode' is deprecated; use `mode' instead报错没关系，可以继续跑
p3 <- plot_cell_trajectory(cds, color_by = "Pseudotime",size=1,show_backbone = TRUE)
p4 <- plot_cell_trajectory(cds, color_by = "celltype",size=1,show_backbone = TRUE)
p5 <- plot_cell_trajectory(cds, color_by = "State",size=1,show_backbone = TRUE)

pdf("./data/output/monocle2_ordergenes.pdf")
p1
p2
p3
p4
p5
dev.off()

time_diff<-differentialGeneTest(cds[ordergene,],cores = 1,fullModelFormulaStr = "~sm.ns(Pseudotime)")
library(magrittr)
library(tidyverse)

time_genes<-time_diff %>% pull(gene_short_name) %>% as.character()
plot_pseudotime_heatmap(cds[time_genes,],num_clusters = 4,show_rownames = T,return_heatmap = T)


time_genes<-top_n(time_diff,n=100,desc(qval)) %>% pull(gene_short_name) %>% as.character() #提取前100个基因
p6 <- plot_pseudotime_heatmap(cds[time_genes,],num_clusters = 4,show_rownames = T,return_heatmap = T)
pdf("./data/output/monocle2_heatmap.pdf")
p6
dev.off()




#版本2：Only_T####
#https://www.jianshu.com/p/5d6fd4561bc0
library(Seurat)
library(ggplot2)
library(cowplot)
library(scater)
library(scran)
library(BiocParallel)
library(BiocNeighbors)
library(data.table)
library(dplyr)
library(Matrix)
library(clustree)
library(monocle)
library(ggsci)
library(ggpubr)

# Only_T<-subset(hms_cluster_id, idents=c("CD4 Naive","CD8 Effect","CD4 Treg",
#                                         "CD8 Exhausted","CD4 Effector Memory","NKT","CD8 Cytotoxic","CD4 Recently Activated","CD8 MAIT",
#                                         "CD8 Progenitor","CD4 Memory","CD4 Effect","CD4 Chronic Activation","CD8 Chronic Activation"))
# DimPlot(Only_T, reduction = "umap", label = TRUE, pt.size = 0.01) 
# DimPlot(Only_T, reduction = "umap", label = FALSE, pt.size = 0.01) 
# saveRDS(Only_T, file = "Only_T_cluster_id_test.rds")


# setwd("~/gse205506/finalize")
# hms_cluster_id<-readRDS("Only_T_cluster_id_test.rds")
# #cd4 cells
# Only_T<-subset(hms_cluster_id, idents=c("CD8 Effect","CD8 Exhausted","CD8 Cytotoxic","CD8 MAIT",
#                                         "CD8 Progenitor","CD8 Chronic Activation"))
# DimPlot(Only_T, reduction = "umap", label = TRUE, pt.size = 0.5) 
# DimPlot(Only_T, reduction = "umap", label = FALSE, pt.size = 0.5) 
# saveRDS(Only_T, file = "cd8_id_test.rds")


setwd("C:/Users/ZFB/Desktop/单细胞生信/GSE212966")

##1.导入seurat对象####
#Only_T
data <- readRDS("./data/temp/Only_T.rds") #（已注释）
# #CD4
# CD4 <- subset(data, idents = c("CD4_Tem", "CD4_Th", "CD4_Tn", "CD4_Treg"))
# data <- CD4
# #CD8
# CD8 <- subset(data, idents = c("CD8_Tc", "CD8_Te", "CD8_Tem", "CD8_Tex", "CD8_Trm"))
# data <- CD8
# #active.ident：细胞类型CD4_Tem
# #tech：Normal/Tumor

#提取表型信息--细胞信息(建议载入细胞的聚类或者细胞类型鉴定信息、实验条件等信息)
expr_matrix <- as(as.matrix(data@assays$RNA@counts), 'sparseMatrix')
#提取表型信息到p_data(phenotype_data)里面 
p_data <- data@meta.data
#整合每个细胞的细胞鉴定信息到p_data里面。如果已经添加则不必重复添加,Normal_ADJ1->CD4_Tem
p_data$celltype <- data@active.ident
#提取基因信息 如生物类型、gc含量等
f_data <-
  data.frame(gene_short_name = row.names(data),
             row.names = row.names(data))
##expr_matrix的行数与f_data的行数相同(gene number), expr_matrix的列数与p_data的行数相同(cell number)

##2.构建CDS对象####
pd <- new('AnnotatedDataFrame', data = p_data) 
fd <- new('AnnotatedDataFrame', data = f_data)

#返回的结果是否为TURE，若不是可以运行下面的代码,否则下一句报错
identical(rownames(fd),rownames(expr_matrix))
expr_matrix<- expr_matrix[rownames(fd), ]
identical(rownames(fd),rownames(expr_matrix))

#将p_data和f_data从data.frame转换AnnotatedDataFrame对象。
cds <- newCellDataSet(expr_matrix,
                      phenoData = pd,
                      featureData = fd,
                      lowerDetectionLimit = 0.5,
                      expressionFamily = negbinomial.size())

#估计size factor和离散度
cds<-estimateSizeFactors(cds)
cds<-estimateDispersions(cds)

##3.过滤低质量的细胞####
cds <- detectGenes(cds,min_expr = 0.1)
print(head(fData(cds)))
#⭐过滤掉在小于10个细胞中表达的基因，原文设置为10个
##设置为0，实际差异不大，相比设置为10，diff基因808 ->801
expressed_genes <- row.names(subset(fData(cds),num_cells_expressed >=10))


#轨迹定义基因选择及可视化和构建轨迹
## 4.轨迹定义基因####
#法1：使用seurat选择的高变基因
express_genes<-VariableFeatures(data)
cds<-setOrderingFilter(cds,express_genes)
plot_ordering_genes(cds)
#法2：使用clusters差异表达基因️⚠️
#deg.cluster<-FindAllMarkers(Epithelial)
#express_genes<-subset(deg.cluster,p_val_adj<0.05)$gene
#cds<-setOrderingFilter(cds,express_genes)
#plot_ordering_genes(cds)
#法3：使用monocle选择的高变基因⚠️
#disp_table<-dispersionTable(cds)
#disp.genes<-subset(disp_table,mean_expression >=0.1 & dispersion_empirical >= 1*dispersion_fit)$gene_id
#cds<-setOrderingFilter(cds,disp.genes)
#plot_ordering_genes(cds)

#~后面是表示对谁做差异分析的变量，理论上可以为p_data的任意列名
diff<-differentialGeneTest(cds[expressed_genes,],fullModelFormulaStr = "~celltype",cores = 1)
head(diff)

##差异表达基因作为轨迹构建的基因,差异基因的选择标准是qval<0.01,decreasing=F表示按数值增加排序
deg<-subset(diff,qval<0.01)
deg<-deg[order(deg$qval,decreasing = F),]
head(deg)
write.table(deg, file = "./data/temp/Only_T_monocle2_deg3.xls",col.names = T,row.names = F,sep = "\t",quote = F)

## 6.轨迹构建基因可视化####
#得到想要的基因列表后，使用setOrderingFilter将它嵌入cds对象，后续的一系列操作都要依赖于这个list。
#setOrderingFilter之后，这些基因被储存在cds@featureData@data[["use_for_ordering"]]，可以通过table(cds@featureData@data[["use_for_ordering"]])查看
ordergene<-rownames(deg)
cds<-setOrderingFilter(cds,ordergene)
cds@featureData@data[["use_for_ordering"]]
table(cds@featureData@data[["use_for_ordering"]])
# FALSE  TRUE 
# 1193   807 
# 保留了807个基因
plot_ordering_genes(cds)
#出的图黑色的点表示用来构建轨迹的差异基因，灰色表示背景基因。
#红色的线是根据第2步计算的基因表达大小和离散度分布的趋势(可以看到，找到的基因属于离散度比较高的基因)

##7.使用DDRTree降維####
#反向图嵌入(DDRTree)算法
cds<-reduceDimension(cds,max_components = 2, method='DDRTree')

##8.⭐开始：细胞排序和轨迹构建####
cds<- orderCells(cds)#时间较久
#⚠️使用root_state参数可以设置拟时间轴的根，如下面的拟时间着色图中可以看出，左边是根。根据state图可以看出，根是State1，若要想把另一端设为根，可以按如下操作
#cds <- orderCells(cds, root_state = 5) #把State5设成拟时间轴的起始点
## 以下报错没关系，可以继续跑
## In dfs(graph = graph, root = root, mode = mode, unreachable = unreachable,  :Argument `neimode' is deprecated; use `mode' instead
saveRDS(cds, file = "./data/temp/Only_T_cds.rds")
cds <- readRDS("./data/temp/Only_T_cds.rds")#不一致

#⚠️m默认排序有点随机的感觉，需要选择和调整
# cds <- orderCells(cds, root_state = 3)#不太对
cds <- orderCells(cds, root_state = 2)#与第一次的结果一致，CD4_Tn在前
# cds <- orderCells(cds, root_state = 1)#不太对，和不设置是一样的
saveRDS(cds, file = "./data/temp/Only_T_cds_root_state=2.rds")
cds <- readRDS("./data/temp/Only_T_cds_root_state=2.rds")

##9.拟时序图####
pdf("./data/output/Only_T_monocle_拟时序图root_state=2.pdf",width = 6,height = 4.3)
# 位置cds /phenoData/data
plot_cell_trajectory(cds, color_by = "Pseudotime",size=1,show_backbone = TRUE) + facet_wrap("~tech",nrow = 1)
plot_cell_trajectory(cds, color_by = "celltype",size=1,show_backbone = TRUE)+ facet_wrap("~tech",nrow = 1)
plot_cell_trajectory(cds, color_by = "State",size=1,show_backbone = TRUE)+ facet_wrap("~tech",nrow = 1)
# p5 <- plot_cell_trajectory(cds, color_by = "cell",size=1,show_backbone = TRUE)
plot_cell_trajectory(cds, color_by = "tech",size=1,show_backbone = TRUE)
#拆分轨迹显示
plot_cell_trajectory(cds, color_by = "State") + facet_wrap("~State",nrow = 1)+ scale_color_npg()
# p8 <- plot_cell_trajectory(cds, color_by = "celltype") + facet_wrap("~cell",nrow = 1)+scale_color_nejm()
# colour=c("#DC143C","#0000FF","#20B2AA","#FFA500","#9370DB","#98FB98","#14a0dc","#b414dc","#123456","#654321")
# p9 <- plot_cell_trajectory(cds, color_by = "celltype") + facet_wrap("~cell",nrow = 1)+scale_color_manual(values = colour)
dev.off()

##10.树形图####
pdf("./data/output/Only_T_monocle_树形图root_state=2.pdf",width = 6,height = 6)
#总体观
plot_cell_trajectory(cds, x=1,y=2,color_by = "celltype") +
  theme(legend.position = 'none',panel.border = element_blank())#+去掉第一個的legend
  #scale_color_manual(values = colour)
plot_complex_cell_trajectory(cds, x=1,y=2,color_by = "celltype") + 
  #scale_color_manual(values = colour)+
  theme(legend.title = element_blank())
#Tumor/Normal对比
plot_cell_trajectory(cds, x=1,y=2,color_by = "celltype") + facet_wrap("~tech",nrow = 1)+
  theme(legend.position = 'none',panel.border = element_blank())#+去掉第一個的legend
  #scale_color_manual(values = colour)
plot_complex_cell_trajectory(cds, x=1,y=2,color_by = "celltype") +facet_wrap("~tech",nrow = 1)+
  #scale_color_manual(values = colour)+
  theme(legend.title = element_blank())
dev.off()

##11.细胞密度图####

###⭐NKT独立####
df<-pData(cds)
# View(df)
###1.ALL####
ClusterName_color_fill<-c("CD4_Tem"="#FBBAB6","CD4_Th"="#EBC77F","CD4_Tn"="#D1D27F", "CD4_Treg"="#9CDA7F","CD8_Tc"="#7FDFBE","CD8_Te"="#7FDFE1","CD8_Tem"="#7FD7FA","CD8_Tex"="#CAC7FF","CD8_Trm"="#F3B5F9","NKT"="#FFB0DD")
ClusterName_color_panel<-c("CD4_Tem"="#F8766D","CD4_Th"="#D89000","CD4_Tn"="#9CC023", "CD4_Treg"="#39B600","CD8_Tc"="#00BF7D","CD8_Te"="#00BFC4","CD8_Tem"="#00B0F6","CD8_Tex"="#9590FF","CD8_Trm"="#E76BF3","NKT"="#FF62BC")
p1 <- ggplot(df, aes(Pseudotime, colour = celltype, fill = celltype)) + facet_wrap("~tech", nrow = 1) +
  geom_density(bw = 0.5, size = 1, alpha = 0.5) + theme_classic2() +
  scale_fill_manual(name = "", values = ClusterName_color_fill) +
  scale_color_manual(name = "", values = ClusterName_color_panel)
###2.CD4_Treg CD4_Tn CD4_Th####
ClusterName_color_fill<-c("CD4_Tem"="transparent","CD4_Th"="#EBC77F","CD4_Tn"="#D1D27F", "CD4_Treg"="#9CDA7F","CD8_Tc"="transparent","CD8_Te"="transparent","CD8_Tem"="transparent","CD8_Tex"="transparent","CD8_Trm"="transparent","NKT"="transparent")
ClusterName_color_panel<-c("CD4_Tem"="transparent","CD4_Th"="#D89000","CD4_Tn"="#9CC023", "CD4_Treg"="#39B600","CD8_Tc"="transparent","CD8_Te"="transparent","CD8_Tem"="transparent","CD8_Tex"="transparent","CD8_Trm"="transparent","NKT"="transparent")
p2 <- ggplot(df, aes(Pseudotime, colour = celltype, fill = celltype)) + facet_wrap("~tech", nrow = 1) +
  geom_density(bw = 0.5, size = 1, alpha = 0.5) + theme_classic2() +
  scale_fill_manual(name = "", values = ClusterName_color_fill) +
  scale_color_manual(name = "", values = ClusterName_color_panel)
###3.CD4_Tem CD8_Tc（毒性）,CD8_Te####
ClusterName_color_fill<-c("CD4_Tem"="#FBBAB6","CD4_Th"="transparent","CD4_Tn"="transparent", "CD4_Treg"="transparent","CD8_Tc"="#7FDFBE","CD8_Te"="#7FDFE1","CD8_Tem"="transparent","CD8_Tex"="transparent","CD8_Trm"="transparent","NKT"="transparent")
ClusterName_color_panel<-c("CD4_Tem"="#F8766D","CD4_Th"="transparent","CD4_Tn"="transparent", "CD4_Treg"="transparent","CD8_Tc"="#00BF7D","CD8_Te"="#00BFC4","CD8_Tem"="transparent","CD8_Tex"="transparent","CD8_Trm"="transparent","NKT"="transparent")
p3 <- ggplot(df, aes(Pseudotime, colour = celltype, fill = celltype)) + facet_wrap("~tech", nrow = 1) +
  geom_density(bw = 0.5, size = 1, alpha = 0.5) + theme_classic2() +
  scale_fill_manual(name = "", values = ClusterName_color_fill) +
  scale_color_manual(name = "", values = ClusterName_color_panel)
###4.CD8_Tem CD8_Tex CD8_Trm####
ClusterName_color_fill<-c("CD4_Tem"="transparent","CD4_Th"="transparent","CD4_Tn"="transparent", "CD4_Treg"="transparent","CD8_Tc"="transparent","CD8_Te"="transparent","CD8_Tem"="#7FD7FA","CD8_Tex"="#CAC7FF","CD8_Trm"="#F3B5F9","NKT"="transparent")
ClusterName_color_panel<-c("CD4_Tem"="transparent","CD4_Th"="transparent","CD4_Tn"="transparent", "CD4_Treg"="transparent","CD8_Tc"="transparent","CD8_Te"="transparent","CD8_Tem"="#00B0F6","CD8_Tex"="#9590FF","CD8_Trm"="#E76BF3","NKT"="transparent")
p4 <- ggplot(df, aes(Pseudotime, colour = celltype, fill = celltype)) + facet_wrap("~tech", nrow = 1) +
  geom_density(bw = 0.5, size = 1, alpha = 0.5) + theme_classic2() +
  scale_fill_manual(name = "", values = ClusterName_color_fill) +
  scale_color_manual(name = "", values = ClusterName_color_panel)
###5.NKT####
ClusterName_color_fill<-c("CD4_Tem"="transparent","CD4_Th"="transparent","CD4_Tn"="transparent", "CD4_Treg"="transparent","CD8_Tc"="transparent","CD8_Te"="transparent","CD8_Tem"="transparent","CD8_Tex"="transparent","CD8_Trm"="transparent","NKT"="#FFB0DD")
ClusterName_color_panel<-c("CD4_Tem"="transparent","CD4_Th"="transparent","CD4_Tn"="transparent", "CD4_Treg"="transparent","CD8_Tc"="transparent","CD8_Te"="transparent","CD8_Tem"="transparent","CD8_Tex"="transparent","CD8_Trm"="transparent","NKT"="#FF62BC")
p5 <- ggplot(df, aes(Pseudotime, colour = celltype, fill = celltype)) + facet_wrap("~tech", nrow = 1) +
  geom_density(bw = 0.5, size = 1, alpha = 0.5) + theme_classic2() +
  scale_fill_manual(name = "", values = ClusterName_color_fill) +
  scale_color_manual(name = "", values = ClusterName_color_panel)

pdf("./data/output/Only_T_monocle2_细胞密度图1root_state=2.pdf",width = 12,height = 2)
p1
p2
p3
p4
p5
dev.off()



###⭐NKT合并   √  ####
df<-pData(cds)
# View(df)
###1.ALL####
ClusterName_color_fill<-c("CD4_Tem"="#FBBAB6","CD4_Th"="#EBC77F","CD4_Tn"="#D1D27F", "CD4_Treg"="#9CDA7F","CD8_Tc"="#7FDFBE","CD8_Te"="#7FDFE1","CD8_Tem"="#7FD7FA","CD8_Tex"="#CAC7FF","CD8_Trm"="#F3B5F9","NKT"="#FFB0DD")
ClusterName_color_panel<-c("CD4_Tem"="#F8766D","CD4_Th"="#D89000","CD4_Tn"="#9CC023", "CD4_Treg"="#39B600","CD8_Tc"="#00BF7D","CD8_Te"="#00BFC4","CD8_Tem"="#00B0F6","CD8_Tex"="#9590FF","CD8_Trm"="#E76BF3","NKT"="#FF62BC")
p1 <- ggplot(df, aes(Pseudotime, colour = celltype, fill = celltype)) + facet_wrap("~tech", nrow = 1) +
  geom_density(bw = 0.5, size = 1, alpha = 0.5) + theme_classic2() +
  scale_fill_manual(name = "", values = ClusterName_color_fill) +
  scale_color_manual(name = "", values = ClusterName_color_panel)
###2.CD4_Treg CD4_Tn CD4_Th####
ClusterName_color_fill<-c("CD4_Tem"="transparent","CD4_Th"="#EBC77F","CD4_Tn"="#D1D27F", "CD4_Treg"="#9CDA7F","CD8_Tc"="transparent","CD8_Te"="transparent","CD8_Tem"="transparent","CD8_Tex"="transparent","CD8_Trm"="transparent","NKT"="transparent")
ClusterName_color_panel<-c("CD4_Tem"="transparent","CD4_Th"="#D89000","CD4_Tn"="#9CC023", "CD4_Treg"="#39B600","CD8_Tc"="transparent","CD8_Te"="transparent","CD8_Tem"="transparent","CD8_Tex"="transparent","CD8_Trm"="transparent","NKT"="transparent")
p2 <- ggplot(df, aes(Pseudotime, colour = celltype, fill = celltype)) + facet_wrap("~tech", nrow = 1) +
  geom_density(bw = 0.5, size = 1, alpha = 0.5) + theme_classic2() +
  scale_fill_manual(name = "", values = ClusterName_color_fill) +
  scale_color_manual(name = "", values = ClusterName_color_panel)
###3.CD4_Tem CD8_Tc（毒性） NKT####
ClusterName_color_fill<-c("CD4_Tem"="#FBBAB6","CD4_Th"="transparent","CD4_Tn"="transparent", "CD4_Treg"="transparent","CD8_Tc"="#7FDFBE","CD8_Te"="transparent","CD8_Tem"="transparent","CD8_Tex"="transparent","CD8_Trm"="transparent","NKT"="#FFB0DD")
ClusterName_color_panel<-c("CD4_Tem"="#F8766D","CD4_Th"="transparent","CD4_Tn"="transparent", "CD4_Treg"="transparent","CD8_Tc"="#00BF7D","CD8_Te"="transparent","CD8_Tem"="transparent","CD8_Tex"="transparent","CD8_Trm"="transparent","NKT"="#FF62BC")
p3 <- ggplot(df, aes(Pseudotime, colour = celltype, fill = celltype)) + facet_wrap("~tech", nrow = 1) +
  geom_density(bw = 0.5, size = 1, alpha = 0.5) + theme_classic2() +
  scale_fill_manual(name = "", values = ClusterName_color_fill) +
  scale_color_manual(name = "", values = ClusterName_color_panel)
###4.CD8_Te CD8_Tem CD8_Tex CD8_Trm####
ClusterName_color_fill<-c("CD4_Tem"="transparent","CD4_Th"="transparent","CD4_Tn"="transparent", "CD4_Treg"="transparent","CD8_Tc"="transparent","CD8_Te"="#7FDFE1","CD8_Tem"="#7FD7FA","CD8_Tex"="#CAC7FF","CD8_Trm"="#F3B5F9","NKT"="transparent")
ClusterName_color_panel<-c("CD4_Tem"="transparent","CD4_Th"="transparent","CD4_Tn"="transparent", "CD4_Treg"="transparent","CD8_Tc"="transparent","CD8_Te"="#00BFC4","CD8_Tem"="#00B0F6","CD8_Tex"="#9590FF","CD8_Trm"="#E76BF3","NKT"="transparent")
p4 <- ggplot(df, aes(Pseudotime, colour = celltype, fill = celltype)) + facet_wrap("~tech", nrow = 1) +
  geom_density(bw = 0.5, size = 1, alpha = 0.5) + theme_classic2() +
  scale_fill_manual(name = "", values = ClusterName_color_fill) +
  scale_color_manual(name = "", values = ClusterName_color_panel)

pdf("./data/output/Only_T_monocle2_细胞密度图2root_state=2.pdf",width = 12,height = 2)
p1
p2
p3
p4
dev.off()






##12.指定基因的可视化化####
###（未使用）基因随细胞状态等的表达变化####
keygenes<-head(ordergene,4)
cds_subset<-cds[keygenes,]
plot_genes_in_pseudotime(cds_subset,color_by = "State")
plot_genes_in_pseudotime(cds_subset,color_by = "celltype")
plot_genes_in_pseudotime(cds_subset,color_by = "Pseudotime")
#指定基因
#CD4_Tem
s.genes <- c("TNFAIP3", "KLRB1","S100A4")#（AQP3的是标志基因OS没问题，但表达太低）
#CD4_Tn
s.genes <- c("CCR7","LTB","GPR183")#👌CCR7
#CD4_Th
s.genes <- c("RGS1", "CTLA4")#（GK 的OS没问题，表达太低而且不是标志基因）
#CD4_Treg
s.genes <- c("BATF", "TIGIT","CTLA4","IL2RA","TNFRSF4")##👌TNFRSF4 

#CD8_Tc
s.genes <- c("CCL5","CCL4","CRTAM","ZNF331")#👌ZNF331
#CD8_Te
s.genes <- c("GZMK", "GZMA","DUSP2")#
#CD8_Tem
s.genes <- c("CRTAM", "DUSP4")#
#CD8_Tex
s.genes <- c("CCL4L2")#👌CCL4L2  #"IFNG",X
#CD8_Trm
s.genes <- c("GZMH")#"NKG7"

#NKT
s.genes <- c("GNLY", "NKG7","KLRD1","PRF1","FGFBP2")#👌FGFBP2

#merge
s.genes <- c("TNFAIP3", "KLRB1",#（AQP3的是标
             "CCR7","LTB","GPR183",#👌CCR7
             "RGS1", "CTLA4",#（GK 的OS没问题，表达太
             "BATF", "TIGIT")
s.genes <- c("CTLA4","TNFRSF4",##👌TNFRSF4 
             "CCL5","CCL4","CRTAM","ZNF331",#👌ZNF331
             "GZMK", "GZMA","DUSP2")#
s.genes <- c("CRTAM", "DUSP4",#
             "CCL4L2",#👌CCL4L2
             "GZMH",#"NKG7"
             "GNLY", "NKG7","KLRD1","PRF1","FGFBP2")#👌FGFBP2
s.genes <- c("EBI3")#"NKG7"
cds[s.genes,]
p1 <- plot_genes_in_pseudotime(cds[s.genes,],color_by = "celltype")
p2 <- plot_genes_jitter(cds[s.genes,],color_by = "celltype")
p3 <- plot_genes_violin(cds[s.genes,],color_by = "celltype")
pdf("./data/output/Only_T_monocle2_指定基因相关表达前last10.pdf",width = 6,height = 8)
plot_genes_in_pseudotime(cds[s.genes,],color_by = "celltype")
# p1
# p2
# p3
dev.off()

###1.基因随细胞状态等的表达变化####
features<- c("TNFAIP3",                               #CD4_Tem
             "CCR7","LTB","GPR183",                   #CD4_Tn  👌CCR7
             "RGS1", "CTLA4",                         #CD4_Th
             "BATF", "TIGIT","IL2RA","TNFRSF4",       #CD4_Treg👌TNFRSF4 
             "CCL5","CCL4","CRTAM","ZNF331",          #CD8_Tc  👌ZNF331
             "GZMK", "GZMA","DUSP2",                  #CD8_Te
             "DUSP4",                                 #CD8_Tem 
             "CCL4L2",                                #CD8_Tex 👌CCL4L2
             "NKG7","GZMH",                           #CD8_Trm   
             "GNLY","KLRD1","PRF1","FGFBP2")          #NKT     👌FGFBP2
plot_genes_in_pseudotime(cds[features,],color_by="celltype",ncol=5)

###3.指定基因绘制热图####
Time_genes <- c("TNFAIP3",                               #CD4_Tem
                "CCR7","LTB","GPR183",                   #CD4_Tn  👌CCR7
                "RGS1", "CTLA4",                         #CD4_Th
                "BATF", "TIGIT","IL2RA","TNFRSF4",       #CD4_Treg👌TNFRSF4 
                "CCL5","CCL4","CRTAM","ZNF331",          #CD8_Tc  👌ZNF331
                "GZMK", "GZMA","DUSP2",                  #CD8_Te
                "DUSP4",                                 #CD8_Tem 
                "CCL4L2",                                #CD8_Tex 👌CCL4L2
                "NKG7","GZMH",                           #CD8_Trm   
                "GNLY","KLRD1","PRF1","FGFBP2")          #NKT     👌FGFBP2
pdf("./data/output/Only_T_monocle2_指定基因相关表达_热图_合并.pdf",width = 6,height = 6)
plot_pseudotime_heatmap(cds[Time_genes,],num_clusters = 3,show_rownames = T,return_heatmap = T)
dev.off()


###2.拟时序展示单个基因表达量####
colnames(pData(cds))
pData(cds)$FGFBP2 =log2(exprs(cds)['FGFBP2',]+1)
plot_cell_trajectory(cds,color_by = "FGFBP2")+scale_color_gsea()

pData(cds)$TNFRSF4 =log2(exprs(cds)['TNFRSF4',]+1)
plot_cell_trajectory(cds,color_by = "TNFRSF4")+scale_color_gsea()

pData(cds)$ZNF331 =log2(exprs(cds)['ZNF331',]+1)
plot_cell_trajectory(cds,color_by = "ZNF331")+scale_color_gsea()


##13.寻找拟时相关的基因####
deg <- read.table(file = "./data/temp/Only_T_monocle2_deg3.xls",header = TRUE)
ordergene <- rownames(deg) 
cds <- setOrderingFilter(cds, ordergene)  
nrow(cds)
nrow(ordergene)
Time_diff<-differentialGeneTest(cds[ordergene,],cores = 1,fullModelFormulaStr = "~sm.ns(Pseudotime)")
Time_diff<-Time_diff[,c(5,2,3,4,1,6,7)]

Time_genes<-Time_diff %>% pull(gene_short_name) %>%  as.character()
plot_pseudotime_heatmap(cds[Time_genes,],num_clusters = 4,show_rownames = T,return_heatmap = T)
p=plot_pseudotime_heatmap(cds[Time_genes,],num_clusters = 4,show_rownames = T,return_heatmap = T)
write.csv(Time_diff,"Time_diff")
#前面通過設置num_clusters將熱圖聚成四個cluster,把每一個cluster的基因單獨提出來分析
p$tree_row
clusters<-cutree(p$tree_row,k=4)
clustering<-data.frame(clusters)
clustering[,1]<-as.character(clustering[,1])
colnames(clustering)<-"Gene_Clusters"
table(clustering)
#提取前100個
Time_genes<-top_n(Time_diff,n=100,desc(qval)) %>% pull(gene_short_name) %>% as.character()
p=plot_pseudotime_heatmap(cds[Time_genes,],num_clusters = 4,show_rownames = T,return_heatmap = T)
#顯著差異基因按熱圖結果排序並保存
hp.genes<-p$tree_row$table[p$tree_row$order]
Time_diff_sig<-Time_diff[hp.genes,c("gene_short_name","pval","qval")]
write.csv(Time_diff_sig,"Time_diff_sig.csv",row.names = F)
p

#提取前50個
Time_genes<-top_n(Time_diff,n=50,desc(qval)) %>% pull(gene_short_name) %>% as.character()
p=plot_pseudotime_heatmap(cds[Time_genes,],num_clusters = 4,show_rownames = T,return_heatmap = T)
#顯著差異基因按熱圖結果排序並保存
hp.genes<-p$tree_row$table[p$tree_row$order]
Time_diff_sig<-Time_diff[hp.genes,c("gene_short_name","pval","qval")]
write.csv(Time_diff_sig,"Time_diff_sig.csv",row.names = F)
p

#提取前23個
Time_genes<-top_n(Time_diff,n=23,desc(qval)) %>% pull(gene_short_name) %>% as.character()
p=plot_pseudotime_heatmap(cds[Time_genes,],num_clusters = 4,show_rownames = T,return_heatmap = T)
#顯著差異基因按熱圖結果排序並保存
hp.genes<-p$tree_row$table[p$tree_row$order]
Time_diff_sig<-Time_diff[hp.genes,c("gene_short_name","pval","qval")]
write.csv(Time_diff_sig,"Time_diff_sig.csv",row.names = F)
p


features<-c("ADH1C","EGR1","SULT1B1","C2CD4A")
features<-c("NAMPT","FABP5","ATP1B1","LGALS4","TUBA1B",
            "ADH1C","EGR1","SULT1B1","ETHE1","EXOC6")
features<-c("HOPX")
features<-c("HSPA1A")

features<-c("TUBA1B")
features<-c("HIST2H2BE")
VlnPlot(,features=features,group.by = "cell")
plot_genes_in_pseudotime(cds[features,],color_by="celltype",ncol=5)




#指定基因
s.genes<-c("ANPEP","CA2")
a<-data.frame(Time_genes)
write.csv(a,"a")

features<-c("HOPX","HSPA1A", "GIMAP4","MT2A","SLC2A3")
VlnPlot(Epithelial,features=features,group.by = "cell")
plot_genes_in_pseudotime(cds[features,],color_by="celltype")
features<-c("HOPX")


#版本2：Duct_epithelial_cell####
#https://www.jianshu.com/p/5d6fd4561bc0
library(Seurat)
library(ggplot2)
library(cowplot)
library(scater)
library(scran)
library(BiocParallel)
library(BiocNeighbors)
library(data.table)
library(dplyr)
library(Matrix)
library(clustree)
library(monocle)
library(ggsci)
library(ggpubr)
setwd("C:/Users/ZFB/Desktop/单细胞生信/GSE212966")

##1.导入seurat对象####
#Duct_epithelial_cell
data <- readRDS("./data/temp/Duct_epithelial_cell_cluster_id_test.rds")
levels(x = data)
data <- readRDS("./data/temp/Only_T.rds") #（已注释）
table(data@meta.data$labels)#singleR之前注释的名字
table(data@meta.data$tech)
table(data@meta.data$celltype)
Idents(data)
levels(data)
#labels存放的自动注释的分类名，这里替换掉，细胞分群名，其他分群也这样修改⭐
data$labels = as.character(Idents(data))
table(data@meta.data$labels) 
saveRDS(data, "./data/temp/Duct_epithelial_cell_cluster_id_test.rds")

#提取表型信息--细胞信息(建议载入细胞的聚类或者细胞类型鉴定信息、实验条件等信息)
expr_matrix <- as(as.matrix(data@assays$RNA@counts), 'sparseMatrix')
#提取表型信息到p_data(phenotype_data)里面 
p_data <- data@meta.data
#整合每个细胞的细胞鉴定信息到p_data里面。如果已经添加则不必重复添加,Normal_ADJ1->CD4_Tem
p_data$celltype <- data@active.ident
#提取基因信息 如生物类型、gc含量等
f_data <-
  data.frame(gene_short_name = row.names(data),
             row.names = row.names(data))
##expr_matrix的行数与f_data的行数相同(gene number), expr_matrix的列数与p_data的行数相同(cell number)

##2.构建CDS对象####
pd <- new('AnnotatedDataFrame', data = p_data) 
fd <- new('AnnotatedDataFrame', data = f_data)

#返回的结果是否为TURE，若不是可以运行下面的代码,否则下一句报错
identical(rownames(fd),rownames(expr_matrix))
expr_matrix<- expr_matrix[rownames(fd), ]
identical(rownames(fd),rownames(expr_matrix))

#将p_data和f_data从data.frame转换AnnotatedDataFrame对象。
cds <- newCellDataSet(expr_matrix,
                      phenoData = pd,
                      featureData = fd,
                      lowerDetectionLimit = 0.5,
                      expressionFamily = negbinomial.size())

#估计size factor和离散度
cds<-estimateSizeFactors(cds)
cds<-estimateDispersions(cds)

##3.过滤低质量的细胞####
cds <- detectGenes(cds,min_expr = 0.1)
print(head(fData(cds)))
#⭐过滤掉在小于10个细胞中表达的基因，原文设置为10个
##设置为0，实际差异不大，相比设置为10，diff基因808 ->801
expressed_genes <- row.names(subset(fData(cds),num_cells_expressed >=10))


#轨迹定义基因选择及可视化和构建轨迹
## 4.轨迹定义基因####
#法1：使用seurat选择的高变基因
express_genes<-VariableFeatures(data)
cds<-setOrderingFilter(cds,express_genes)
plot_ordering_genes(cds)
#法2：使用clusters差异表达基因️⚠️
#deg.cluster<-FindAllMarkers(Epithelial)
#express_genes<-subset(deg.cluster,p_val_adj<0.05)$gene
#cds<-setOrderingFilter(cds,express_genes)
#plot_ordering_genes(cds)
#法3：使用monocle选择的高变基因⚠️
#disp_table<-dispersionTable(cds)
#disp.genes<-subset(disp_table,mean_expression >=0.1 & dispersion_empirical >= 1*dispersion_fit)$gene_id
#cds<-setOrderingFilter(cds,disp.genes)
#plot_ordering_genes(cds)

Track_genes<- graph_test(cds,neighbor_graph="principal_graph",cores=6)
Track_genes <-Track_genes[,c(5,2,3,4,1,6)]%>%filter(q_value <1e-3)
write.csv(Track_genes,"Trajectory_genes.csv", row.names = F)


#~后面是表示对谁做差异分析的变量，理论上可以为p_data的任意列名
diff<-differentialGeneTest(cds[expressed_genes,],fullModelFormulaStr = "~celltype",cores = 1)
head(diff)

##差异表达基因作为轨迹构建的基因,差异基因的选择标准是qval<0.01,decreasing=F表示按数值增加排序
deg<-subset(diff,qval<0.01)
deg<-deg[order(deg$qval,decreasing = F),]
head(deg)
write.table(deg, file = "./data/temp/Duct_epithelial_cell_monocle2_deg3.xls",col.names = T,row.names = F,sep = "\t",quote = F)

## 6.轨迹构建基因可视化####
#得到想要的基因列表后，使用setOrderingFilter将它嵌入cds对象，后续的一系列操作都要依赖于这个list。
#setOrderingFilter之后，这些基因被储存在cds@featureData@data[["use_for_ordering"]]，可以通过table(cds@featureData@data[["use_for_ordering"]])查看
ordergene<-rownames(deg)
cds<-setOrderingFilter(cds,ordergene)
cds@featureData@data[["use_for_ordering"]]
table(cds@featureData@data[["use_for_ordering"]])
# FALSE  TRUE 
# 220  1780 
# 保留了1780个基因
plot_ordering_genes(cds)
#出的图黑色的点表示用来构建轨迹的差异基因，灰色表示背景基因。
#红色的线是根据第2步计算的基因表达大小和离散度分布的趋势(可以看到，找到的基因属于离散度比较高的基因)

##7.使用DDRTree降維####
#反向图嵌入(DDRTree)算法
cds<-reduceDimension(cds,max_components = 2, method='DDRTree')

##8.⭐开始：细胞排序和轨迹构建####
cds<- orderCells(cds)#时间较久
#⚠️使用root_state参数可以设置拟时间轴的根，如下面的拟时间着色图中可以看出，左边是根。根据state图可以看出，根是State1，若要想把另一端设为根，可以按如下操作
#cds <- orderCells(cds, root_state = 5) #把State5设成拟时间轴的起始点
## 以下报错没关系，可以继续跑
## In dfs(graph = graph, root = root, mode = mode, unreachable = unreachable,  :Argument `neimode' is deprecated; use `mode' instead
saveRDS(cds, file = "./data/temp/Duct_epithelial_cell_cds.rds")
cds <- readRDS("./data/temp/Duct_epithelial_cell_cds.rds")#不一致

#⚠️m默认排序有点随机的感觉，需要选择和调整
# cds <- orderCells(cds, root_state = 3)#不太对
cds <- orderCells(cds, root_state = 1)#与第一次的结果一致，CD4_Tn在前
# cds <- orderCells(cds, root_state = 1)#不太对，和不设置是一样的
saveRDS(cds, file = "./data/temp/Duct_epithelial_cell_cds_root_state=1.rds")
cds <- readRDS("./data/temp/Duct_epithelial_cell_cds_root_state=1.rds")

##9.拟时序图####
pdf("./data/output/Duct_epithelial_cell_monocle_拟时序图root_state=1.pdf",width = 6,height = 4.3)
# 位置cds /phenoData/data
plot_cell_trajectory(cds, color_by = "Pseudotime",size=1,show_backbone = TRUE) + facet_wrap("~tech",nrow = 1)
plot_cell_trajectory(cds, color_by = "celltype",size=1,show_backbone = TRUE)+ facet_wrap("~tech",nrow = 1)
plot_cell_trajectory(cds, color_by = "State",size=1,show_backbone = TRUE)+ facet_wrap("~tech",nrow = 1)
# p5 <- plot_cell_trajectory(cds, color_by = "cell",size=1,show_backbone = TRUE)
plot_cell_trajectory(cds, color_by = "tech",size=1,show_backbone = TRUE)
#拆分轨迹显示
plot_cell_trajectory(cds, color_by = "State") + facet_wrap("~State",nrow = 1)+ scale_color_npg()
# p8 <- plot_cell_trajectory(cds, color_by = "celltype") + facet_wrap("~cell",nrow = 1)+scale_color_nejm()
# colour=c("#DC143C","#0000FF","#20B2AA","#FFA500","#9370DB","#98FB98","#14a0dc","#b414dc","#123456","#654321")
# p9 <- plot_cell_trajectory(cds, color_by = "celltype") + facet_wrap("~cell",nrow = 1)+scale_color_manual(values = colour)
dev.off()

##10.树形图####
pdf("./data/output/Duct_epithelial_cell_monocle_树形图root_state=1.pdf",width = 6,height = 6)
#总体观
plot_cell_trajectory(cds, x=1,y=2,color_by = "celltype") +
  theme(legend.position = 'none',panel.border = element_blank())#+去掉第一個的legend
#scale_color_manual(values = colour)
plot_complex_cell_trajectory(cds, x=1,y=2,color_by = "celltype") + 
  #scale_color_manual(values = colour)+
  theme(legend.title = element_blank())
#Tumor/Normal对比
plot_cell_trajectory(cds, x=1,y=2,color_by = "celltype") + facet_wrap("~tech",nrow = 1)+
  theme(legend.position = 'none',panel.border = element_blank())#+去掉第一個的legend
#scale_color_manual(values = colour)
plot_complex_cell_trajectory(cds, x=1,y=2,color_by = "celltype") +facet_wrap("~tech",nrow = 1)+
  #scale_color_manual(values = colour)+
  theme(legend.title = element_blank())
dev.off()

##11.细胞密度图####

df<-pData(cds)
# View(df)
###1.ALL####
ClusterName_color_fill<-c("Group_1"="#FBBAB6","Group_2"="#EBC77F","Group_3"="#D1D27F", "Group_4"="#9CDA7F","Group_5"="#7FDFBE","Group_6"="#7FDFE1","Group_7"="#7FD7FA","Group_8"="#CAC7FF","Group_9"="#F3B5F9")
ClusterName_color_panel<-c("Group_1"="#F8766D","Group_2"="#D89000","Group_3"="#9CC023", "Group_4"="#39B600","Group_5"="#00BF7D","Group_6"="#00BFC4","Group_7"="#00B0F6","Group_8"="#9590FF","Group_9"="#E76BF3")
p1 <- ggplot(df, aes(Pseudotime, colour = celltype, fill = celltype)) + facet_wrap("~tech", nrow = 1) +
  geom_density(bw = 0.5, size = 1, alpha = 0.5) + theme_classic2() +
  scale_fill_manual(name = "", values = ClusterName_color_fill) +
  scale_color_manual(name = "", values = ClusterName_color_panel)
###2.Group_1249####
ClusterName_color_fill<-c("Group_1"="#FBBAB6","Group_2"="#EBC77F","Group_3"="transparent", "Group_4"="#9CDA7F","Group_5"="transparent","Group_6"="transparent","Group_7"="transparent","Group_8"="transparent","Group_9"="#F3B5F9")
ClusterName_color_panel<-c("Group_1"="#F8766D","Group_2"="#D89000","Group_3"="transparent", "Group_4"="#39B600","Group_5"="transparent","Group_6"="transparent","Group_7"="transparent","Group_8"="transparent","Group_9"="#E76BF3")
p2 <- ggplot(df, aes(Pseudotime, colour = celltype, fill = celltype)) + facet_wrap("~tech", nrow = 1) +
  geom_density(bw = 0.5, size = 1, alpha = 0.5) + theme_classic2() +
  scale_fill_manual(name = "", values = ClusterName_color_fill) +
  scale_color_manual(name = "", values = ClusterName_color_panel)
###3.Group_35678####
ClusterName_color_fill<-c("Group_1"="transparent","Group_2"="transparent","Group_3"="#D1D27F", "Group_4"="transparent","Group_5"="#7FDFBE","Group_6"="#7FDFE1","Group_7"="#7FD7FA","Group_8"="#CAC7FF","Group_9"="transparent")
ClusterName_color_panel<-c("Group_1"="transparent","Group_2"="transparent","Group_3"="#9CC023", "Group_4"="transparent","Group_5"="#00BF7D","Group_6"="#00BFC4","Group_7"="#00B0F6","Group_8"="#9590FF","Group_9"="transparent")
p3 <- ggplot(df, aes(Pseudotime, colour = celltype, fill = celltype)) + facet_wrap("~tech", nrow = 1) +
  geom_density(bw = 0.5, size = 1, alpha = 0.5) + theme_classic2() +
  scale_fill_manual(name = "", values = ClusterName_color_fill) +
  scale_color_manual(name = "", values = ClusterName_color_panel)


pdf("./data/output/Duct_epithelial_cell_monocle2_细胞密度图1root_state=1.pdf",width = 12,height = 2)
p1
p2
p3
dev.off()

##12.指定基因的可视化化####
###（未使用）基因随细胞状态等的表达变化####
keygenes<-head(ordergene,4)
cds_subset<-cds[keygenes,]
plot_genes_in_pseudotime(cds_subset,color_by = "State")
plot_genes_in_pseudotime(cds_subset,color_by = "celltype")
plot_genes_in_pseudotime(cds_subset,color_by = "Pseudotime")
#指定基因
#CD4_Tem
s.genes <- c("TNFAIP3", "KLRB1","S100A4")#（AQP3的是标志基因OS没问题，但表达太低）
#CD4_Tn
s.genes <- c("CCR7","LTB","GPR183")#👌CCR7
#CD4_Th
s.genes <- c("RGS1", "CTLA4")#（GK 的OS没问题，表达太低而且不是标志基因）
#CD4_Treg
s.genes <- c("BATF", "TIGIT","CTLA4","IL2RA","TNFRSF4")##👌TNFRSF4 

#CD8_Tc
s.genes <- c("CCL5","CCL4","CRTAM","ZNF331")#👌ZNF331
#CD8_Te
s.genes <- c("GZMK", "GZMA","DUSP2")#
#CD8_Tem
s.genes <- c("CRTAM", "DUSP4")#
#CD8_Tex
s.genes <- c("CCL4L2")#👌CCL4L2  #"IFNG",X
#CD8_Trm
s.genes <- c("GZMH")#"NKG7"

#NKT
s.genes <- c("GNLY", "NKG7","KLRD1","PRF1","FGFBP2")#👌FGFBP2

#merge
s.genes <- c("TNFAIP3", "KLRB1",#（AQP3的是标
             "CCR7","LTB","GPR183",#👌CCR7
             "RGS1", "CTLA4",#（GK 的OS没问题，表达太
             "BATF", "TIGIT")
s.genes <- c("CTLA4","TNFRSF4",##👌TNFRSF4 
             "CCL5","CCL4","CRTAM","ZNF331",#👌ZNF331
             "GZMK", "GZMA","DUSP2")#
s.genes <- c("CRTAM", "DUSP4",#
             "CCL4L2",#👌CCL4L2
             "GZMH",#"NKG7"
             "GNLY", "NKG7","KLRD1","PRF1","FGFBP2")#👌FGFBP2
s.genes <- c("EBI3")#"NKG7"
cds[s.genes,]
p1 <- plot_genes_in_pseudotime(cds[s.genes,],color_by = "celltype")
p2 <- plot_genes_jitter(cds[s.genes,],color_by = "celltype")
p3 <- plot_genes_violin(cds[s.genes,],color_by = "celltype")
pdf("./data/output/Only_T_monocle2_指定基因相关表达前last10.pdf",width = 6,height = 8)
plot_genes_in_pseudotime(cds[s.genes,],color_by = "celltype")
# p1
# p2
# p3
dev.off()

###1.基因随细胞状态等的表达变化####
features<- c("TNFAIP3",                               #CD4_Tem
             "CCR7","LTB","GPR183",                   #CD4_Tn  👌CCR7
             "RGS1", "CTLA4",                         #CD4_Th
             "BATF", "TIGIT","IL2RA","TNFRSF4",       #CD4_Treg👌TNFRSF4 
             "CCL5","CCL4","CRTAM","ZNF331",          #CD8_Tc  👌ZNF331
             "GZMK", "GZMA","DUSP2",                  #CD8_Te
             "DUSP4",                                 #CD8_Tem 
             "CCL4L2",                                #CD8_Tex 👌CCL4L2
             "NKG7","GZMH",                           #CD8_Trm   
             "GNLY","KLRD1","PRF1","FGFBP2")          #NKT     👌FGFBP2
plot_genes_in_pseudotime(cds[features,],color_by="celltype",ncol=5)

###3.指定基因绘制热图####
#早
Time_1 <- c("CEACAM7", "MUC3A", "HS3ST1", "GABRP", "LINC01133", "GALNT5", "CLDN18", 
            "MUCL3", "CORO2A", "ITGA2", "MGLL", "MISP",
            "RAB11FIP1", "ITGA6", "NDRG1", "LDLR", "KLF5", "CEACAM6", "KLK10",
            "PLAT", "CEACAM5", "KCNK1", "SQLE", "CENPF", "HMGA1", "EMP2", "CDKN2A",
            "PLAAT3", "MALL", "KLF4", "PLAC8", "MAL2")         
#中
Time_2 <- c("TCN1", "CRISP3", "PGC", "SERPINA1", "MUC5B", "PIGR", "AREG",
            "TFF2", "TPM2", "AHNAK2", "TFF3", "MMP1", "GPX2", "IGFBP2", "TFF1", "LYZ", "CYP3A5", "GCNT3", "ARL14", "DMBT1", "ANXA10",
            "MLPH", "CXCL17", "ITGB6", "DUOX2", "CAPN8", "MECOM", "SMIM24", "ZG16B", "TMC5", "SLC44A4", "SFN",
            "KRT17", "CTSE", "CCL20", "TSPAN1", "TGM2", "AGR2", "IFI27", "LCN2", "LAMC2", "MUC1", "MARCKSL1", "LGALS3",
            "JPT1", "C19orf33", "PHLDA2", "MSLN", "GPRC5A", "ELF3", "CAMK2N1", "NQO1", "LGALS4",
            "LAMB3", "TMPRSS4", "S100P", "FOXQ1", "APOL1", "CXCL5",
            "FXYD3", "S100A4", "KRT7", "SLPI", "CSTB", "KRT19", "S100A6", "CRIP1", "C15orf48", "PLAUR", "CLDN4", "S100A11")          #NKT     👌FGFBP2
#晚
Time_3 <- c("REG1A", "GP2", "REG1B", "PNLIPRP1", "PRSS2", "SCTR",
            "KLK1", "CTRC", "CELA3B", "CEL", "SYCN", "CPA2", "CELA2A", "MUC6", "CTRB1",
            "PNLIP", "CTRB2", "PLA2G1B", "CPB1", "CPA1", "CELA3A", "CLPS", "PRSS1", "CLU",
            "TTYH1", "VTCN1", "RBP1", "MT1G", "DEFB1", "LEFTY1", "ITIH5", "HEG1", "TPM1",
           "IL1R1", "AQP1", "SERPINA5", "RASSF4", "SERPINA6", "CFTR", "SLC4A4", "FGFR3", "C6", "SLC3A1", 
            "TTN", "GUCY1A1", "FXYD2", "CITED4", "FILIP1L",
            "AMBP", "COL18A1", "CCND1", "SH3YL1", "ENC1", "GATM", "CHST9", "CLDN10", "GMNN", "DCDC2", "SORBS2",
            "EPHB6", "MEG3", "ID4", "HES4", "FGFR2", "PROX1", "SFRP5", "SOD3", "CADM1"
          )          #NKT     👌FGFBP2

plot_pseudotime_heatmap(cds[Time_3,], show_rownames = T,num_clusters = 1)

pdf("./data/output/Duct_epithelial_cell_monocle2_heatmap_Time1.pdf",width = 12)
plot_pseudotime_heatmap(cds[Time_1,], show_rownames = T,num_clusters = 1)#num_clusters = 3
dev.off()
pdf("./data/output/Duct_epithelial_cell_monocle2_heatmap_Time2.pdf",width = 12)
plot_pseudotime_heatmap(cds[Time_2,], show_rownames = T,num_clusters = 1)
dev.off()
pdf("./data/output/Duct_epithelial_cell_monocle2_heatmap_Time3.pdf",width = 12)
plot_pseudotime_heatmap(cds[Time_3,], show_rownames = T,num_clusters = 1)
dev.off()


###2.拟时序展示单个基因表达量####
colnames(pData(cds))
pData(cds)$FGFBP2 =log2(exprs(cds)['FGFBP2',]+1)
plot_cell_trajectory(cds,color_by = "FGFBP2")+scale_color_gsea()

pData(cds)$TNFRSF4 =log2(exprs(cds)['TNFRSF4',]+1)
plot_cell_trajectory(cds,color_by = "TNFRSF4")+scale_color_gsea()

pData(cds)$ZNF331 =log2(exprs(cds)['ZNF331',]+1)
plot_cell_trajectory(cds,color_by = "ZNF331")+scale_color_gsea()



#13####
cds <- readRDS("./data/temp/Duct_epithelial_cell_cds_root_state=1.rds")
print(head(fData(cds)))
diff<-differentialGeneTest(cds[expressed_genes,],fullModelFormulaStr = "~celltype",cores = 1)
deg<-subset(diff,qval<0.01)
deg<-deg[order(deg$qval,decreasing = F),]
write.table(deg, file = "./data/temp/Duct_epithelial_cell_monocle2_deg.xls",col.names = T,row.names = F,sep = "\t",quote = F)

ordergene<-rownames(deg)
cds<-setOrderingFilter(cds,ordergene)
cds<-reduceDimension(cds,max_components = 2, method='DDRTree')

time_diff<-differentialGeneTest(cds[ordergene,],cores = 1,fullModelFormulaStr = "~sm.ns(Pseudotime)")
library(magrittr)
library(tidyverse)

# time_genes<-time_diff %>% pull(gene_short_name) %>% as.character()
# plot_pseudotime_heatmap(cds[time_genes,],num_clusters = 4,show_rownames = T,return_heatmap = T)
time_genes_100<-top_n(time_diff,n=100,desc(qval)) %>% pull(gene_short_name) %>% as.character() #提取前100个基因
time_genes_200<-top_n(time_diff,n=200,desc(qval)) %>% pull(gene_short_name) %>% as.character() #提取前100个基因
time_genes_300<-top_n(time_diff,n=300,desc(qval)) %>% pull(gene_short_name) %>% as.character() #提取前100个基因

deg_top_100<-time_diff[time_genes_100,]
deg_top_200<-time_diff[time_genes_200,]
deg_top_300<-time_diff[time_genes_300,]

deg_top_100<-deg_top_100[order(deg_top_100$num_cells_expressed,decreasing = T),]
deg_top_200<-deg_top_200[order(deg_top_200$num_cells_expressed,decreasing = T),]
deg_top_300<-deg_top_300[order(deg_top_300$num_cells_expressed,decreasing = T),]

write.csv(deg_top_100,file = "./data/output/Duct_epithelial_cell_monocle2_deg100.csv")
write.csv(deg_top_200,file = "./data/output/Duct_epithelial_cell_monocle2_deg200.csv")
write.csv(deg_top_300,file = "./data/output/Duct_epithelial_cell_monocle2_deg300.csv")

p1 <- plot_pseudotime_heatmap(cds[time_genes_100,],num_clusters = 2,show_rownames = T,return_heatmap = T)
p2 <- plot_pseudotime_heatmap(cds[time_genes_100,],num_clusters = 3,show_rownames = T,return_heatmap = T)
p3 <- plot_pseudotime_heatmap(cds[time_genes_100,],num_clusters = 4,show_rownames = T,return_heatmap = T)

pdf("./data/output/Duct_epithelial_cell_monocle2_heatmap_deg100-2.pdf",width = 12,height = 12)
p1
dev.off()
pdf("./data/output/Duct_epithelial_cell_monocle2_heatmap_deg100-3.pdf",width = 12,height = 12)
p2
dev.off()
pdf("./data/output/Duct_epithelial_cell_monocle2_heatmap_deg100-4.pdf",width = 12,height = 12)
p3
dev.off()



saveRDS(test, file = "test_monocle.rds")
write.table(pseudotime_de, file = "pseudotime_de.rds", quote = FALSE, sep = '\t', row.names = FALSE, col.names = TRUE)
write.table(states_de, file = "states_de.rds", quote = FALSE, sep = '\t', row.names = FALSE, col.names = TRUE)


#13.寻找拟时相关的基因####
cds <- readRDS("./data/temp/Duct_epithelial_cell_cds_root_state=1.rds")
deg <- read.table(file = "./data/temp/Duct_epithelial_cell_monocle2_deg3.xls",header = TRUE)
ordergene <- rownames(deg) 
cds <- setOrderingFilter(cds, ordergene)  
nrow(cds)
nrow(ordergene)
Time_diff<-differentialGeneTest(cds[ordergene,],cores = 1,fullModelFormulaStr = "~sm.ns(Pseudotime)")
Time_diff<-differentialGeneTest(cds,cores = 1,fullModelFormulaStr = "~sm.ns(Pseudotime)")

Time_diff<-Time_diff[,c(5,2,3,4,1,6,7)]

Time_genes<-Time_diff %>% pull(gene_short_name) %>%  as.character()
plot_pseudotime_heatmap(cds[Time_genes,],num_clusters = 3,show_rownames = F,return_heatmap = T)
p=plot_pseudotime_heatmap(cds[Time_genes,],num_clusters = 4,show_rownames = T,return_heatmap = T)
write.csv(Time_diff,"Time_diff")
#前面通過設置num_clusters將熱圖聚成四個cluster,把每一個cluster的基因單獨提出來分析
p$tree_row
clusters<-cutree(p$tree_row,k=4)
clustering<-data.frame(clusters)
clustering[,1]<-as.character(clustering[,1])
colnames(clustering)<-"Gene_Clusters"
table(clustering)
#提取前100個
Time_genes<-top_n(Time_diff,n=100,desc(qval)) %>% pull(gene_short_name) %>% as.character()
p=plot_pseudotime_heatmap(cds[Time_genes,],num_clusters = 4,show_rownames = T,return_heatmap = T)
#顯著差異基因按熱圖結果排序並保存
hp.genes<-p$tree_row$table[p$tree_row$order]
Time_diff_sig<-Time_diff[hp.genes,c("gene_short_name","pval","qval")]
write.csv(Time_diff_sig,"Time_diff_sig.csv",row.names = F)
p

#提取前50個
Time_genes<-top_n(Time_diff,n=50,desc(qval)) %>% pull(gene_short_name) %>% as.character()
p=plot_pseudotime_heatmap(cds[Time_genes,],num_clusters = 4,show_rownames = T,return_heatmap = T)
#顯著差異基因按熱圖結果排序並保存
hp.genes<-p$tree_row$table[p$tree_row$order]
Time_diff_sig<-Time_diff[hp.genes,c("gene_short_name","pval","qval")]
write.csv(Time_diff_sig,"Time_diff_sig.csv",row.names = F)
p

#提取前23個
Time_genes<-top_n(Time_diff,n=23,desc(qval)) %>% pull(gene_short_name) %>% as.character()
p=plot_pseudotime_heatmap(cds[Time_genes,],num_clusters = 4,show_rownames = T,return_heatmap = T)
#顯著差異基因按熱圖結果排序並保存
hp.genes<-p$tree_row$table[p$tree_row$order]
Time_diff_sig<-Time_diff[hp.genes,c("gene_short_name","pval","qval")]
write.csv(Time_diff_sig,"Time_diff_sig.csv",row.names = F)
p


features<-c("ADH1C","EGR1","SULT1B1","C2CD4A")
features<-c("NAMPT","FABP5","ATP1B1","LGALS4","TUBA1B",
            "ADH1C","EGR1","SULT1B1","ETHE1","EXOC6")
features<-c("HOPX")
features<-c("HSPA1A")

features<-c("TUBA1B")
features<-c("HIST2H2BE")
VlnPlot(,features=features,group.by = "cell")
plot_genes_in_pseudotime(cds[features,],color_by="celltype",ncol=5)




#指定基因
s.genes<-c("ANPEP","CA2")
a<-data.frame(Time_genes)
write.csv(a,"a")

features<-c("HOPX","HSPA1A", "GIMAP4","MT2A","SLC2A3")
VlnPlot(Epithelial,features=features,group.by = "cell")
plot_genes_in_pseudotime(cds[features,],color_by="celltype")
features<-c("HOPX")